#include "CentralColumnSolver.h"
#include <iostream>
#include <omp.h>
#include <algorithm>
#include "Perfs.h"

using namespace std;
CentralColumnSolver::CentralColumnSolver()
{
    //ctor
}

CentralColumnSolver::~CentralColumnSolver()
{
    //dtor
}

inline
void CentralColumnSolver::update_result(Result & result, t_int x0, t_int y0, t_int x1, t_int y1, t_int sum)
{
    result.x0 = x0;
    result.y0 = y0;
    result.x1 = x1;
    result.y1 = y1;
    result.sum = sum;
}

inline
void CentralColumnSolver::initialize_accumulators( Result & best_result,
                                                  vector<int> & accumulator_left, vector<int> & accumulator_right,
                                                  int left_x, int local_y1, int left_limit_x, int right_limit_x)
{
    const Mat & mat = *mat_ptr;
    const vector <int> & row = mat[local_y1];
    accumulator_left[0] = row[left_x];
    accumulator_right[0] = row[left_x + 1];
    int i = 1;
    //initializing upper structures
    for (int local_x_left = left_x - 1;local_x_left >= left_limit_x ; local_x_left--)
    {
        accumulator_left[i] = accumulator_left[i-1] + row[local_x_left];
        i++;
    }
    //initializing lower structures
    i = 1;
    for (int local_x_right = left_x + 2 ; local_x_right <= right_limit_x; local_x_right++)
    {
        accumulator_right[i] = accumulator_right[i-1] + row[local_x_right];
        i++;
    }

    //testing against actual best result
    int max_left, dist_left, max_right, dist_right;
    find_max_and_index(accumulator_left, max_left, dist_left);
    find_max_and_index(accumulator_right, max_right, dist_right);
    int total_sum = max_left + max_right;
    if (total_sum > best_result.sum)
    {
        update_result(best_result, left_x - dist_left, local_y1,
                           left_x + 1 + dist_right, local_y1, total_sum);
    }

}
inline
void CentralColumnSolver::find_max_and_index(const vector <int > & accumulator, int & maxi,
                                             int & dist_from_center)
{
    const vector<int>::const_iterator max_elem = max_element(accumulator.begin(), accumulator.end());
    maxi = *max_elem;
    dist_from_center = distance(accumulator.begin(), max_elem);
}


void CentralColumnSolver::solve_MCS(const Data * original_data,  Result & result, const s_rect & rect,
                                    int central_column)
{
    this->data_ptr = original_data;
    this->mat_ptr = &(original_data->getMat());
    const Mat & mat = *mat_ptr;
/*
	cout << "Central col : (" << rect.x0 << "," << rect.y0 << ") to (";
	cout << rect.x1 << "," 	<< rect.y1 << ") with central col = ";
	cout << central_column << "and thread = " << omp_get_thread_num() << endl;
*/
	int nb_cols_left = central_column - rect.x0 + 1;
	int nb_cols_right = rect.x1 - central_column;
	const int left_x = rect.x0 + nb_cols_left - 1;

	//best_result.display();
#pragma omp parallel
{
	Result best_result;

	best_result.sum = 1 <<31;

	vector<int> accumulator_left(nb_cols_left);
	vector<int> accumulator_right(nb_cols_right);
    #pragma omp for schedule(static, 10)
	for (int local_y1 = rect.y0 ; local_y1 <= rect.y1 ; local_y1++)
	{
        //cout << "go" << endl;
	    initialize_accumulators( best_result,
                             accumulator_left, accumulator_right,
                              left_x, local_y1,rect.x0, rect.x1);
		const int limit_y = rect.y1;
		for (int local_y2 = local_y1 + 1 ; local_y2 <= limit_y ; local_y2++)
		{
           int local_accumulator_left = 0;
           int i = 0;
           const int limit_left_x = rect.x0;
		   for (int local_x_left = left_x ; local_x_left >= limit_left_x ; local_x_left-- )
		   {
		       //cout << "x = " << local_x_left << " : " << accumulator_left[i] << endl;
				local_accumulator_left += mat[local_y2][local_x_left];
				accumulator_left[i] += local_accumulator_left;
               // cout << "x = " << local_x_left << " : " << accumulator_left[i] << endl;
				i++;
			}
           i = 0;
        int local_accumulator_right = 0;
        const int limit_right_x = rect.x1;
		   for (int local_x_right = left_x + 1 ; local_x_right <= limit_right_x ; local_x_right++ )
		   {
				local_accumulator_right += mat[local_y2][local_x_right];
				accumulator_right[i] += local_accumulator_right;
				i++;
			}
			int max_left, dist_left_from_center, max_right, dist_right_from_center;
			find_max_and_index(accumulator_left, max_left, dist_left_from_center);
			find_max_and_index(accumulator_right, max_right, dist_right_from_center);
			int total_sum = max_left + max_right;
			if (total_sum > best_result.sum)
			{
			    //cout << "new best result found" << endl;
			    //best_result.update(left_result.x0, left_result.y0, right_result.x1, right_result.y1, total_sum);
			    best_result.x0 = left_x - dist_left_from_center;
			    best_result.y0 = local_y1;
			    best_result.x1 = left_x + 1 + dist_right_from_center;
			    best_result.y1 = local_y2;
			    best_result.sum = total_sum;
			    //best_result.display();
			}
		}
	}
	#pragma omp critical
	{
        if (best_result > result)
            result = best_result;
	}
    //best_result.display();
}//pragma omp parallel
   // if (best_result > result)
     //   result = best_result;
    //best_result.display();
    //final_result.display();
    //result = final_result;
}

inline
void CentralColumnSolver::initialize_matrix(int local_y1, int left_x, const s_rect & rect, vector< vector <int> > & accumulator_left,
                  vector< vector <int> > &accumulator_right, Result & best_result)
{
    const Mat & mat = *(this->mat_ptr);
    #pragma omp parallel for schedule(static,10)
for (int local_y2 = local_y1 ; local_y2 <= rect.y1 ; local_y2++)
		{
			int local_accumulator_left = 0;
           int i = 0;
		   for (int local_x_left = left_x ; local_x_left >= rect.x0 ; local_x_left-- )
		   {
		       //cout << "x = " << local_x_left << " : " << accumulator_left[i] << endl;
				local_accumulator_left += mat[local_y2][local_x_left];
				accumulator_left[local_y2][i] += local_accumulator_left;
				i++;
			}

           i = 0;
           //initialise best lower result with x1, left_x+1, x2, left_x+1
            int local_accumulator_right = 0;
		   for (int local_x_right = left_x + 1 ; local_x_right <= rect.x1 ; local_x_right++ )
		   {
				local_accumulator_right += mat[local_y2][local_x_right];
				accumulator_right[local_y2][i] += local_accumulator_right;
				i++;
			}

			int max_left, max_right, dist_left, dist_right;
			find_max_and_index(accumulator_left[local_y2], max_left, dist_left);
			find_max_and_index(accumulator_right[local_y2], max_right, dist_right);
			int total_sum = max_left + max_right;
			if (total_sum > best_result.sum)
			{
			    update_result(best_result, left_x - dist_left, local_y1,
                     left_x + 1 + dist_right, local_y2, total_sum);
			}
		}
}
inline
void CentralColumnSolver::substract_row_with_row( vector <int> & plus_row, const vector<int> & minus_row)
{
    const int v_size = plus_row.size();
    for (int column = 0; column < v_size ; column++)
    {
        plus_row[column] -= minus_row[column];
    }
}
inline
void CentralColumnSolver::update_matrix(int local_y1, int left_x, const s_rect & rect,
                                         vector< vector <int> > & accumulator_left,
                  vector< vector <int> > &accumulator_right, Result & best_result)
{
    int first_row = local_y1+1;
    #pragma omp parallel for schedule(static,10)
    for(int row = local_y1 +1; row < rect.y1 ; row++)
    {
        int acc_row = row - local_y1;
        const vector<int> & origin_left_row = accumulator_left[local_y1];
        const vector<int> & origin_right_row = accumulator_right[local_y1];
        vector<int> & update_left_row = accumulator_left[acc_row];
        vector<int> & update_right_row = accumulator_right[acc_row];
        substract_row_with_row(update_left_row, origin_left_row);
        substract_row_with_row(update_right_row, origin_right_row);

			int max_left, max_right, dist_left, dist_right;
			find_max_and_index(update_left_row, max_left, dist_left);
			find_max_and_index(update_right_row, max_right, dist_right);
        int total_sum = max_left + max_right;
        if (total_sum > best_result.sum)
        {
            int actual_left_col = left_x - dist_left;
            int actual_right_col = left_x + 1 + dist_right;
            update_result(best_result, actual_left_col, first_row, actual_right_col, row, total_sum);
        }

    }
}

void CentralColumnSolver::solve_MCS_optim(const Data * original_data,  Result & result, const s_rect & rect,
                                    int central_column)
{
    this->data_ptr = original_data;
    this->mat_ptr = &(original_data->getMat());
    const Mat & mat = *mat_ptr;
    //#warning "Test what happens when trying to solve on a single col"
/*
	cout << "Central col : (" << rect.x0 << "," << rect.y0 << ") to (";
	cout << rect.x1 << "," 	<< rect.y1 << ") with central col = ";
	cout << central_column << "and thread = " << omp_get_thread_num() << endl;
*/
    int nb_rows = rect.y1 - rect.y0 + 1;
	int nb_cols_left = central_column - rect.x0 + 1;
	int nb_cols_right = rect.x1 - central_column;
	int left_x = rect.x0 + nb_cols_left - 1;
    result.sum = 1 << 31;
	//best_result.display();

   // #pragma omp parallel
    {

double time1 = Perfs::my_gettimeofday();
	Result best_result;

	update_result(best_result,left_x, rect.y0, left_x+1, rect.x0,
                    mat[rect.y0][left_x] + mat[rect.y0][left_x +1]);

	vector < vector<int> > accumulator_left(nb_rows);
    vector < vector<int> > accumulator_right(nb_rows);
	for (int row = 0 ; row < nb_rows ; row++)
        accumulator_left[row].resize(nb_cols_left);
    for (int row = 0 ; row < nb_rows ; row++)
        accumulator_right[row].resize(nb_cols_right);
//#pragma omp barrier
double time2 = Perfs::my_gettimeofday();
cout << time2-time1 << " to allocate mem" << endl;
    //#pragma omp for schedule(static, 10)
	for (int local_y1 = rect.y0 ; local_y1 <= rect.y1 ; local_y1++)
	{
	    //initialise matrix of sums from local_y1 to rect.y1
	    //double init1 = Perfs::my_gettimeofday();
	    initialize_matrix(local_y1, left_x, rect, accumulator_left, accumulator_right
                       , best_result);
        //double mid = Perfs::my_gettimeofday();
        update_matrix(local_y1, left_x, rect, accumulator_left, accumulator_right
                       , best_result);
        //double endi = Perfs::my_gettimeofday();
        //cout << "int = " << mid-init1 << " and " << endi-mid << endl;

    }
    //#pragma omp critical
    {
        if (best_result > result)
            result = best_result;
    }
       }//omp parallel

}

